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ABSTRACT 

Light transmission through circular subwavelength apertures in metallic films with surrounding nanostructures 
is investigated numerically. Numerical results are obtained with a frequency-domain finite-clement method. 
Convergence of the obtained observables to very low levels of numerical error is demonstrated. Very good 
agreement to experimental results from the literature is reached, and the utility of the method is demonstrated 
in the investigation of the influence of geometrical parameters on enhanced transmission through the apertures. 
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1. INTRODUCTION 

Experiments by Lezec et al. investigating light transmission through circular subwavelength apertures in metallic 
films with surrounding nanostructures have revealed surprising effects of enhanced transmission and a collimated 
beam of transmitted light. 1 The mechanisms involved in the transmission spectra can be understood using 2D 
models and approximative methods. 2:3 However, for a fully quantitative understanding and from the point of 
view of optics design, accurate numerical simulations of Maxwell's equations for the 3D, cylindrically symmetric 
problem are desired. Popov et al. have developed a well converging field expansion method with a Fourier-Bessel 
basis in order to simulate similar setups. 4-6 To our knowledge results on a simulation of the experiment of Lezec 
et al. with this method have not been published. Baida et al. have investigated the setting numerically using a 
finite-difference time-domain method (FDTD), however, an agreement of simulation results with the experiment 
is not reached with this method. 7 

We report on a dedicated finite-element method for the simulation task at hand. With this, accurately con- 
verged results can be obtained in relatively short computation times. We discuss simulations of the experimental 
setup 1 and present simulation results agreeing very well with results of the experiment. 

The paper is structured as follows: Section |2~T1 reviews the general field of light transmission through periodic 
and isolated subwavelength apertures and motivates the need for accurate simulation tools. The investigated 
setup of a single aperture surrounded by circular grooves is defined in Section 12.21 Section 12.31 introduces the 
finite-element solver JCMsuite used in this study. Numerical results are presented in Chapter [3] including a 
convergence study for the investigated setup in Section [3HJ results on the experimental setup of Lezec et al. and 
on the setup investigated by Baida et al. in Section [3. 21 . and results on the influence of a geometrical parameter 
on transmission spectra through the subwavelength apertures in Section [3.31 Chapter [4] concludes the paper. 
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2. BACKGROUND 

2.1 Enhanced transmission through isolated and periodic subwavelength apertures 

In 1997 experiments were reported 8 claiming surprisingly high transmission of light through 2D-periodic sub- 
wavelength hole arrays drilled in a metallic sheet of silver. The enhancement was reported to be orders of 
magnitude larger as compared with the transmission of one single hole of same size predicted by standard 
aperture theory. 9 The explanation of the experiments was mainly based on the existence of surface plasmons 
polaritons (SPPs) that are excited by the incident light field. The experiment has created a lot of attention, and 
the perception of surface plasmons has been widely supported by further investigation. However, in 2004 Lezec 
and Thio report on the experimental observation of both enhancement and suppression in the transmission 
spectra of hole arrays in metallic films. 10 Their work suggests that a lateral interference phenomenon of an 
evanescent wave with the incident light is the driving force in shaping the transmission spectra. Especially, a 
suppression in the transmission spectrum is not be explained by the SPP model which predicts only enhancement 
with respect to the single-hole case. Enhanced transmission with a similar signature is also reported for hole 
arrays in tungsten which does not support SPPs in the investigated regime. The authors further report on 
a transmission enhancement factor not larger than 7, in contrast to the often-quoted enhancement factor of 
order 1000. 10:11 The explanation of the transmission spectra with interference phenomena has been supported 
also by a study showing that the transmission peak positions are determined by the period of the hole array 
and are material independent. 12 Flammer et al. studied transmission through a single slit surrounded by linear 
grating structures. 13 Also this work confirms the explanation through an interference phenomenon. However, 
in this case it has been found that the surface wave is an SPP, and not an evanescent wave. The theoretical 
model of this work shows that a SPP is launched along the metal surface, while interference of the SPP with the 
incident light along with resonant cavity effects give rise to suppression and enhancement in transmission. In 
another work it has even been shown that there are no holes necessary for enhanced transmission through thin 
metallic films provided the film is modulated. 14 

This story of different interpretations of experimental results points out that approximate models often 
trigger in-depth physical understanding of the ongoing physics. However, when different models lead to partly 
contradictive conclusions, a comparison with rigorous simulations has to be performed. These methods need to 
consider the full wavelength dependent material dispersion behaviour, must be able to reach a sufficient spatial 
resolution of the structures and should have a good convergence rate to get fast and accurate results. 

2.2 Investigated Setup 

The investigated setup is relatively simple: A plane wave under normal incidence illuminates a silver sheet which 
is perforated by a small circular hole, surrounded by assist features. The geometry is symmetric under rotation 
around the axis of the hole. Figure [T] schematically shows the geometrical setup, 1 the parameters are given in 
Table [TJ For the relative electric permittivity of silver, e r , we assume a Drude model, e r = 1 — oj 2 / (lj 2 + i oj u> c ), 
with plasma frequency uj p and collision frequency oj c as given in Table [IJ and frequency uj — cq2"k/\q. The 
geometry and material parameters correspond to the parameters of Ref. 7 (structure CAG), investigating the 
same experimental setup. 1 




Figure 1. Schematics of the experimental setup: silver sheet (light blue) of thickness t is perforated by a cylindrical hole of 
diameter d, surrounded by four circular grooves of width w, depth h and distance p. The setup is rotationally symmetric 
around the y-axis. An illuminating plane wave is incident in y-direction. 
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In sets of successive simulations either numerical parameters are scanned in order to investigate the conver- 
gence of the obtained discrete solutions, or physical parameters like the wavelength or geometry parameters are 
scanned in order to obtain, e.g., the spectral response of the structure. The angular spectrum of transmitted light 
is investigated by evaluating Fourier transforms of the computed field at the boundaries of the computational 
domain. For an analysis of the field distributions, these can also be exported and visualized. Figure [2] shows an 
example where a simulated near field is visualized in different 2D sections through the 3D field distribution. The 
sections clearly resemble the cylindrical geometry of the structure. 





set 1 


set 2 


set 3 


set 4 


set 5 


Geometry parameters 


d 


330 nm 


300 nm 


t 
h 


300 nm 
60 nm 


P 


600 nm 


550 nm - 650 nm 


565 nm 


w 


p/2 


Incident spectrum 


A 


670 nm 


400 nm - HOOnm 


500 nm - 800 nm 


631.4 nm 


Material properties 


UJp 

UJ C 


1.374el6 s- 1 
3.21el3 s" 1 



Table 1. Five different sets of used parameter settings, for geometry parameters compare Fig. [TJ Please note that we 
investigate settings with two different values for the central hole diameter d because we compare to results from a numerical 
study 7 using d — 330 nm and we compare to experimental results 1 where the setting d = 300 nm was reported. 
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Figure 2. Pseudocolor images of the magnitude of a simulated electric field, 0\E\, in five cross-sections parallel to the 
surface of the silver film, at different heights y. From left to right: y = —151 nm (1 nm below the surface on the entrance 
side), y = — 91 nm (lnm below the bottom of the grooves on the entrance side), y = Onm (at the center of the film), 
y — 91 nm (1 nm above the bottom of the grooves on the exit side), y = 151 nm (1 nm above the surface on the exit side). 
The pseudocolor images (red: high magnitude, blue: low magnitude) are scaled by a factor j3 relative to the amplitude of 
the incoming plane wave (from left to right: /3 — 4%, 8%, 8%, 20%, 10%). The parameters of the simulation are given in 
Table [T] (set 5). [See original paper for better image resolution.] 

2.3 Numerical Method 

For solving the time-harmonic Maxwell's equations describing light scattering we use a dedicated finite-element 
method for nanooptics simulations: The finite-element package JCMsuite is a joint development of the Zusc 
Institute Berlin and JCMwave. This solver has been successfully applied to a wide range of electromagnetic 
field computations including metamaterials, 15 photonic crystal fibers, 16 nearficld- microscopy, 17 and optical mi- 
crolithography. 18 The solver is also used for pattern reconstruction in optical metrology. 19 The performance of 
our methods has been pointed out in several benchmarks. 20, 21 The main ingredients for an accurate performance 
are adaptive methods, based on goal-oriented error estimators, higher-order vector elements, and fast direct and 
iterative numerical methods for solving matrix equations. 22 

In this contribution we apply the solver for light scattering computations on rotationally symmetric geome- 
tries, included in the programme package JCMsuite. Assuming harmonic time-dependence of the electric field 
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with frequency u>, E t = E exp(— iwt) , Maxwell's equations can be written as 

curl /x" 1 curl E - lu 2 eE = , (la) 

dive.E = . (lb) 

The spatial distribution of the tensorial permittivity e and permeability fi can be chosen (nearly) arbitrarily in 
the computational window Cl- mi C R 3 . The scattering problem is defined as follows: Given an incoming electric 
field E lnc compute the total electric field E (satisfying ([I]) in R 3 ) such that the scattered field E sc = E — E- mc 
defined on the exterior domain f2 oxt = R 3 \ f2i n t is purely outward radiating. 

Using the transformation rules for differential forms with Jacobian matrix J of the transformation (r, y, (f>) <— > 
(x,y,z) and curl* = (d r , d y , d^Y x, E* = 3 t E, £* = |J|J _1 eJ~ t and /i* = |J|J~VJ~S Maxwell's equations in 
cylinder coordinates (r, y, ip) read as 

curl*^ 1 curl — ui 2 e*E. t = , (2a) 
div *e*E* = . (2b) 

We assume now that the material distribution is cylindrically symmetric, i.e., e* = e*(r, y), /i* = /i*(r, y), and 
we expand the transformed incoming field -Ei nc .* = J* -Einc in a Fourier series, 

n— oo 

E inc ,*(r,y,<p)= e mc,«(r,?/)e" lv . 

n— — oo 

With this, the 3D Maxwell's equations separate into 2D problems for each n G Z, 

curl*^^ 1 curU !n e„ - oj 2 e*e n = , (3) 

with curl*.„ = (drjdyjin)* 1 x . Equation ((3|) in its weak formulation is solved several times for increasing \n\. 
The computation is stopped automatically when the incoming modes ej nCi „ become sufficiently small. 

3. NUMERICAL RESULTS 

3.1 Convergence 

We demonstrate sound convergence properties of our finite-element solver by performing several computations of 
the same physical setup with increasing numerical resolution. The physical parameters for this convergence study 
are given in Table [T] (parameter set 1), the wavelength of the incident plane wave corresponds to a transmission 
peak for this parameter setting. 

The physical observable of interest in this case is the far field intensity transmitted through the aperture, T, 
detected at a specific detection angle. T is obtained from the FEM near field solution. The far field intensity 
T(N) computed from a numerical approximation with N degrees of freedom, En, in general differs from T 
computed from the exact solution. The relative error of T(N) is defined as AT(N) = \T(N) ~ T qo |/|T qe |, 
where T qG is the intensity computed from the quasi-exact solution, i.e., from a solution on a finer mesh than the 
meshes of the solutions corresponding to T(N). Indeed, it would be desirable to compare T(N) to an analytical 
solution. However, for problems where an analytical solution is not available, the quasi-exact solution is used as 
a makeshift. To validate the results we have further tested that results of our method converge towards results 
obtained from Mie-Theory (for scattering off a sphere), and we have further quantitatively reproduced published 
results on light transmission through circular holes obtained from a Fourier-Bessel modal method. 4 

Figure [3] shows the dependence of the relative error Ar(A^) on the number of degrees of freedom in the 
numerical approximation, N. The detection angle is in normal direction to the silver film in this case. Increasing 
the finite-element polynomial degree and refining the mesh, both leads to an increase in N. And an increase 
in N leads - when the convergent regime is reached - to a decrease of the relative error of the numerical 
approximation. For the displayed results we have used finite-elements of second polynomial order, and we have 
used up to ten successive adaptive mesh refinement steps. Figure G2 further displays the computational effort in 
terms of computation time on a standard workstation. The software has multi-processing capabilities, however, 
for a better comparability we have used a single core of a cpu only for these computations. The convergence plot 
shows that for achieving a relative error of, e.g., one percent a moderate number of unknowns of about TV = 10 5 
is needed, corresponding to computation times well below one minute on a standard personal computer (PC). 
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Figure 3. a) Convergence plot: Relative error of computed transmission through a nanoaperture in dependence on number 
of unknowns N of the FEM problem, b) Computational effort: Computation time in seconds (single-thread computation 
on a standard workstation) in dependence on N (corresponding to the computations displayed in a). 



3.2 Enhanced transmission and beaming effect 

In order to numerically investigate the experimental setting of Lezec et al. we have performed a wavelength scan 
and recorded far-field transmission into different spatial directions. The physical parameters in this case were 
as defined in Table [T] (parameter set 3). These parameters correspond to the published data by Lezec et al. 
Figure H] shows the simulated spectra. Clearly, enhanced transmission around a peak wavelength of A = 660 nm 
is found for low observation angles. This is in very good agreement with the experimental findings. 1 




Figure 4. a) Intensity spectrum of transmitted light for different detection angles, b) Angular intensity spectrum of 
transmitted light at maximum transmission wavelength (A = 660 nm) . 



For reproducing the beaming effect observed in the experiment we further compute the far-field pattern of 
transmitted light at wavelength of maximum transmission through the nanoaperture. Figure |4]b) shows the 
corresponding angular spectrum of diffracted light: Transmitted light emerges from the aperture as a beam with 
small angular divergence of approximately ±3.5 degree (FWHM). Again these results agree very well with the 
experimental findings. 

For comparison to results obtained with a FDTD method we also include an intensity spectrum for the 
geometrical values defined in Table [1] (parameter set 2). The geometry and material parameters correspond to 
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Figure 5. Intensity spectrum of transmitted light at a detection angle of degree. Geometry as defined in parameter set 
1 of Table [JJ Transmission peak around a wavelength of A = 670 nm. 



the parameters of Ref. 7 Figure [5] shows the computed spectrum. We note that this spectrum is significantly 
different from the spectrum computed by Baida et al. using FDTD (Fig. 2b in Ref. 7 ). The spectrum of Ref. 7 
is also in discrepancy with the experimental results of Lezec et al. In Ref. 7 the discrepancy between simulation 
results and experimental findings is attributed to either the difference between the dispersive properties of the 
real metal used in the experiment and the data used in the simulations, or to a possible difference between 
geometrical parameters used in the modeling and those of the real experiment. However, our results with the 
same geometry and material parameters, as shown in Figure [5J agree well with the experimental findings. A 
slight shift (w 10 nm) in transmission peak wavelength is explained by a difference in aperture diameter between 
numerical parameters (d — 330 nm) and experimental setting (d — 300 nm) . This suggests that the discrepancy 
in Ref. 7 is rather caused by insufficient convergence of the used numerical methods. 

3.3 Dependence of enhanced transmission on geometrical parameters 
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Figure 6. Dependence of the transmission spectrum on the spacing p of the circular grooves surrounding the nanoaperture. 
Intensity spectra of transmitted light at a detection angle of degree (a), resp. 4 degree (b). The groove spacing is varied 
between 550 nm (black, peak at A ~ 618.5 nm) and 650 nm (red, peak at A ~ 700.5 nm). 



We have performed a set of numerical experiments where we simulate transmission spectra for different 
spacings p of the grating structure surrounding the aperture (cf. Fig. [I}. All other parameters of the setup 
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remain fixed to set 4 of Table [TJ As numerical parameters, second order finite-elements and two successive 
adaptive refinement steps of the coarse grids were chosen. According to Section 13. f I this setting corresponds to a 
relative error of the transmission values of about one percent. Figure [5] shows how the transmission maxima for 
detection angles degree and 4 degree are shifted to the red side of the spectrum with increasing spacing p. In 
the investigated regime, the spectral position of the maxima depends nearly linearly on spacing p. The significant 
shift of the transmission peak demonstrates the strong influence of the structures surrounding the nanoaperture, 
emphasizing the importance of interference effects for enhanced transmission through subwavelength apertures. 

A maximum of enhanced transmission is observed at a spacing of p — 565 nm and a wavelength of A = 
631.4 nm. Figure [5] shows the near field amplitude in different cross sections through the 3D field distribution 
for these parameters (parameter set 5 of Table [TJ . It is interesting to see how well the field distribution at the 
exit side of the structure resembles the field distribution on the input side, with a relative amplitude scaling of 
(only) about 2.5. 

4. CONCLUSIONS 

We have presented a finite-element method for accurate and fast, rigorous simulations of light propagation in 
cylindrically symmetric, metallic nano-structures. Simulation results of light transmission through structured 
nano-apertures show an excellent agreement with experimental results from the literature. The method has been 
used to investigate the strong dependency of transmission spectra on a geometrical parameter. 

With the presented finite-element solver, a method is at hand which allows for more detailed studies of light 
transmission through subwavelength apertures, and for an efficient and reliable design of metallic nanostructures 
for specific applications, e.g., in sensing or for near field microscopy. 
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